#############################################
##### Replication Code
##### Mobilizing for Prisoner Release During a Pandemic via Perspective-Taking Exercises
##### Mackenzie Israel-Trummel
#############################################



library(multilevel)
library(texreg)
library(tidyverse)
library(reshape2)
library(effects)
library(stargazer)

# Load Data ---------------------------------------------------------------

load(file = "covid.rdata")


##how many people consented? 
table(covid$Q84) ##2045
table(covid$condition) ##2001 went past the demographic/quota block to get randomized into a treatment, 44 screened out for full quotas



##subset to those who took the survey (dropping people not assigned to a treatment)
table(covid$condition)

covid <- subset(covid, condition != "")

names(covid)
table(covid$Q84)


###################
##Clean Moderator Variables for Heterogeneous Treatment Effects: party and carceral contact
###################

###Party ID
table(covid$Q33)
table(covid$Q94.1)
table(covid$Q95.1)

covid$pid7[covid$Q33 == "Democrat" & covid$Q94.1 == "Strong"] <- 0
covid$pid7[covid$Q33 == "Democrat" & covid$Q94.1 == "Not very strong"] <- 1

covid$pid7[covid$Q95.1 == "Closer to Democratic"] <- 2
covid$pid7[covid$Q95.1 == "Don't know"] <- 3
covid$pid7[covid$Q95.1 == "Closer to Republican"] <- 4

covid$pid7[covid$Q33 == "Republican" & covid$Q94.1 == "Not very strong"] <- 5
covid$pid7[covid$Q33 == "Republican" & covid$Q94.1 == "Strong"] <- 6

table(covid$Q33, covid$pid7)
table(covid$Q94.1, covid$pid7)
table(covid$Q95.1, covid$pid7)

###carceral contact

table(covid$Q52.1)
covid$servedtime <- 0
covid$servedtime[covid$Q52.1 == "Yes"] <-1
table(covid$Q52.1, covid$servedtime)

table(covid$Q77)

table(covid$Q83)


covid$proxcontact_close <- 1
covid$proxcontact_close[covid$Q77 == ""] <- 0
covid$proxcontact_close[covid$Q77 == "Yes, an acquaintance"] <- 0
covid$proxcontact_close[covid$Q77 == "No"] <- 0
table(covid$Q77, covid$proxcontact_close)

covid$pers_or_prox <- 0
covid$pers_or_prox[covid$servedtime == 1] <- 1
covid$pers_or_prox[covid$proxcontact_close == 1] <- 1
table(covid$servedtime, covid$pers_or_prox)
table(covid$proxcontact_close, covid$pers_or_prox)
table(covid$pers_or_prox)


#####################
##Clean Dependent Variables
#####################

##DV attitudes

covid$release_all[covid$Q66_1 == "Strongly disagree"] <- 0
covid$release_all[covid$Q66_1 == "Somewhat disagree"] <- 1
covid$release_all[covid$Q66_1 == "Neither agree nor disagree"] <- 2
covid$release_all[covid$Q66_1 == "Somewhat agree"] <- 3
covid$release_all[covid$Q66_1 == "Strongly agree"] <- 4
table(covid$Q66_1, covid$release_all)


covid$release_old[covid$Q66_2 == "Strongly disagree"] <- 0
covid$release_old[covid$Q66_2 == "Somewhat disagree"] <- 1
covid$release_old[covid$Q66_2 == "Neither agree nor disagree"] <- 2
covid$release_old[covid$Q66_2 == "Somewhat agree"] <- 3
covid$release_old[covid$Q66_2 == "Strongly agree"] <- 4
table(covid$Q66_2, covid$release_old)

covid$release_immuno[covid$Q66_3 == "Strongly disagree"] <- 0
covid$release_immuno[covid$Q66_3 == "Somewhat disagree"] <- 1
covid$release_immuno[covid$Q66_3 == "Neither agree nor disagree"] <- 2
covid$release_immuno[covid$Q66_3 == "Somewhat agree"] <- 3
covid$release_immuno[covid$Q66_3 == "Strongly agree"] <- 4
table(covid$Q66_3, covid$release_immuno)

covid$release_nonviolent[covid$Q66_4 == "Strongly disagree"] <- 0
covid$release_nonviolent[covid$Q66_4 == "Somewhat disagree"] <- 1
covid$release_nonviolent[covid$Q66_4 == "Neither agree nor disagree"] <- 2
covid$release_nonviolent[covid$Q66_4 == "Somewhat agree"] <- 3
covid$release_nonviolent[covid$Q66_4 == "Strongly agree"] <- 4
table(covid$Q66_4, covid$release_nonviolent)

covid$release_bail[covid$Q66_5 == "Strongly disagree"] <- 0
covid$release_bail[covid$Q66_5 == "Somewhat disagree"] <- 1
covid$release_bail[covid$Q66_5 == "Neither agree nor disagree"] <- 2
covid$release_bail[covid$Q66_5 == "Somewhat agree"] <- 3
covid$release_bail[covid$Q66_5 == "Strongly agree"] <- 4
table(covid$Q66_5, covid$release_bail)



##create release attitudes index

covid$release <- covid$release_all + covid$release_old + covid$release_immuno + covid$release_nonviolent + covid$release_bail
table(covid$release)

prop.table(table(covid$release))



### Creating Variable for Semibehavioral Outcome

table(covid$Disagree, covid$Reconciled) ##Disagree shows where any coders disagreed. Reconciled are the codes assigned during reconciliation process for those disagreements
covid$Reconciled[covid$Disagree == "Agree" & covid$Mackenzie == 1] <- 1
covid$Reconciled[covid$Disagree == "Agree" & covid$Mackenzie == 2] <- 2
covid$Reconciled[covid$Disagree == "Agree" & covid$Mackenzie == 3] <- 3
covid$Reconciled[covid$Disagree == "Agree" & covid$Mackenzie == 4] <- 4
covid$Reconciled[covid$Disagree == "Agree" & covid$Mackenzie == 5] <- 5
covid$Reconciled[covid$Disagree == "Agree" & covid$Mackenzie == 99] <- 99
table(covid$Reconciled)
table(covid$Disagree, covid$Reconciled)

##create three level semi-behavioral for the analysis

covid$supportrelease_Rec3 <- 0
covid$supportrelease_Rec3[covid$Reconciled == 1] <- 1
covid$supportrelease_Rec3[covid$Reconciled == 2] <- -1 
table(covid$Reconciled, covid$supportrelease_Rec3)
table(covid$supportrelease_Rec3)


##code variable for whether R identifies as White or not (includes those who ID as White and another group)

table(covid$Q25)
covid$white <- 0
covid$white[covid$Q25 == "White"] <- 1
covid$white[covid$Q25 == "White,American Indian or Alaska Native"] <- 1
covid$white[covid$Q25 == "White,American Indian or Alaska Native,Another group"] <- 1
covid$white[covid$Q25 == "White,Asian American"] <- 1
covid$white[covid$Q25 == "White,Black/African American"] <- 1
covid$white[covid$Q25 == "White,Black/African American,American Indian or Alaska Native"] <- 1
covid$white[covid$Q25 == "White,Black/African American,Hispanic/Latino/Latinx"] <- 1
covid$white[covid$Q25 == "White,Black/African American,Hispanic/Latino/Latinx,American Indian or Alaska Native"] <- 1
covid$white[covid$Q25 == "White,Hispanic/Latino/Latinx"] <- 1
covid$white[covid$Q25 == "White,Hispanic/Latino/Latinx,American Indian or Alaska Native"] <- 1
covid$white[covid$Q25 == "White,Hispanic/Latino/Latinx,Asian American"] <- 1
covid$white[covid$Q25 == "White,Hispanic/Latino/Latinx,Asian American,American Indian or Alaska Native,Another group"] <- 1
covid$white[covid$Q25 == "White,Middle Eastern or North African"] <- 1
table(covid$Q25, covid$white)

##Subset Data to Rs who ID as white

covid.white <- subset(covid, white == 1)
table(covid.white$white) ##1535 White respondents

##Set notreat as reference level for conditions
covid.white$condition <- as.factor(covid.white$condition)
covid.white$condition <- relevel(covid.white$condition, ref = "notreat")


##Descriptive Statistics on DVs
summary(covid.white$release_all)
summary(covid.white$release_old)
summary(covid.white$release_immuno)
summary(covid.white$release_nonviolent)
summary(covid.white$release_bail)

prop.table(table(covid.white$release))

table(covid.white$supportrelease_Rec)
prop.table(table(covid.white$supportrelease_Rec))


##Cronbach's alpha
release <- cbind(covid.white$release_all, covid.white$release_old, covid.white$release_immuno, covid.white$release_nonviolent, covid.white$release_bail)
cronbach(release) 
##ALPHA =  0.8908887



##Alphas across conditions:
## subsets by treatment
control <- subset(covid.white, condition == "control")
notreat <- subset(covid.white, condition == "notreat")
t1 <- subset(covid.white, condition == "t1")
t2 <- subset(covid.white, condition == "t2")

release_control <- cbind(control$release_all, control$release_old, control$release_immuno, control$release_nonviolent, control$release_bail)
cronbach(release_control)

release_notreat <- cbind(notreat$release_all, notreat$release_old, notreat$release_immuno, notreat$release_nonviolent, notreat$release_bail)
cronbach(release_notreat)

release_t1 <- cbind(t1$release_all, t1$release_old, t1$release_immuno, t1$release_nonviolent, t1$release_bail)
cronbach(release_t1)

release_t2 <- cbind(t2$release_all, t2$release_old, t2$release_immuno, t2$release_nonviolent, t2$release_bail)
cronbach(release_t2)



###################
## Create Figure 1
###################


id.w <- covid.white$ResponseId
all.w <- covid.white$release_all
old.w <- covid.white$release_old
immuno.w <- covid.white$release_immuno
nonviolent.w <- covid.white$release_nonviolent
bail.w <- covid.white$release_bail


data <- data.frame(id.w, all.w, old.w, immuno.w, nonviolent.w, bail.w)

releasedata <- melt(data, id.vars = c("id.w"), measure.vars = c("all.w", "old.w", "immuno.w", "nonviolent.w", "bail.w"))

releasedata$group[releasedata$variable == "all.w"] = "Release \nall"
releasedata$group[releasedata$variable == "old.w"] = "Elderly"
releasedata$group[releasedata$variable == "immuno.w"] = "Immunocom-\npromised"
releasedata$group[releasedata$variable == "nonviolent.w"] = "Nonviolent"
releasedata$group[releasedata$variable == "bail.w"] = "Awaiting \nbail"

names(releasedata)


aggregate(releasedata$value, by = list(Category = releasedata$variable), FUN = mean, na.rm = TRUE)


#mean plot
pdf("fig1a.pdf", height= 5, width = 7)
ggplot(releasedata, aes(x=group, y=value)) + 
  stat_summary(fun=mean, geom="bar", width = .25, position=position_dodge(1)) + 
  stat_summary(fun.data = mean_cl_normal,  
               geom = "errorbar", width = .15) + 
  theme_bw() +   
  labs(x = "", y = "Mean") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 18), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))
dev.off()


pdf("fig1b.pdf", height= 5, width = 7)
ggplot(covid.white, aes(x=release)) + 
  geom_bar(mapping = aes(x = release, y = ..prop.., group = 1), stat = "count") + 
  labs(x="Support for Release",y="")+
  theme_bw() +
  ylim(0, .15) +
  theme(axis.text=element_text(size=18), 
        axis.title=element_text(size=18)) 

dev.off()



###############
## Treatment Effects: Figure 2
###############


##treatment conditions


aggregate(covid.white$release, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE) ##means by treatment

aggregate(covid.white$release, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(covid.white$release, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[1, 2] ##diff from no treat

aggregate(covid.white$release, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(covid.white$release, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2, 2] ##diff from control

sd(covid.white$release, na.rm = TRUE)

##treatment effect size in standard deviations: sd is 5.80021
##diff between control and t1: 0.1624966 standard deviations
####diff between control and t2: 0.2827467 standard deviations

##diff between no treatment and t1: 0.1451148 standard deviations
####diff between no treatment  and t2: 0.2653648 standard deviations


#t-tests between treatments, notreat as baseline

##Naive t-tests
t.a <- t.test(control$release, notreat$release)
t.b <- t.test(t1$release, notreat$release)
t.c <- t.test(t2$release, notreat$release) 

t.test(t1$release, control$release)
t.test(t2$release, control$release)
t.test(t2$release, t1$release)


##Recover adjusted p values using Holm correction method
pairwise.t.test(covid.white$release, covid.white$condition, p.adjust.method="holm")


##plot effect of  treatment

treateffect.release <- c(0, t.a$estimate[1]- t.a$estimate[2], t.b$estimate[1]- t.b$estimate[2], t.c$estimate[1]- t.c$estimate[2])

lower.release <- c(0, t.a$conf.int[1], t.b$conf.int[1], t.c$conf.int[1])

upper.release <- c(0, t.a$conf.int[2], t.b$conf.int[2], t.c$conf.int[2])


treatment <- c("No\n treatment", "Informational\n Control",  "Perspective:\n egocentric", "Perspective:\n surrogate")

data.release <- as.data.frame(treateffect.release)
data.release$lower <- lower.release
data.release$upper <- upper.release
data.release$treatment <- treatment
data.release$treatment <- factor(data.release$treatment, levels = c("Perspective:\n surrogate", "Perspective:\n egocentric",  "Informational\n Control", "No\n treatment"))
data.release$labels <- round(data.release$treateffect.release, 3)
data.release$labels[data.release$labels == 0] <- NA
#data.release$labels2 <- c("" , "-0.101"," 0.842", "")

p1 <- ggplot(data = data.release, aes(x = treateffect.release, y = treatment)) +
  geom_text(aes(label = labels), nudge_y = 0.2, nudge_x = -0.002, size=6) +
  theme_bw() +
  geom_point(size = 4) +
  geom_errorbarh(aes(y = treatment, xmin = lower, xmax = upper, height=.1),size=.5) +
  xlab("Effect on Support for Releasing Prisoners")  + ylab("") +
  theme(axis.text = element_text(size = 18), #Increase size of axis text
        axis.title = element_text(size = 18))	+  geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +
  #ggtitle(title)+
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, 0), "cm")) +
  annotate("text", x = 0.4, y = 4.15, parse = TRUE, label = "bar(x)==7.869", col = "black", size = 6) 


pdf("fig2a.pdf", height= 5, width = 7)
p1
dev.off()



###############
##semibehavioral outcome: writing to sheriff
###############

aggregate(covid.white$supportrelease_Rec3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)

aggregate(covid.white$supportrelease_Rec3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2]- aggregate(covid.white$supportrelease_Rec3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[1,2] ##diff from no treatment

aggregate(covid.white$supportrelease_Rec3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2]- aggregate(covid.white$supportrelease_Rec3, by = list(Category = covid.white$condition), FUN = mean, na.rm = TRUE)[2,2] ##diff from control


sd(covid.white$supportrelease_Rec3, na.rm = TRUE) #0.614861
##standard deviation changes
##control to t1: 0.144766
##control to t2: 0.1175669
##no treatment to t1: 0.2415746
##no treatment to t2: 0.2143754

#t-tests between treatments, notreat as baseline

##Naive t tests

tb.a <- t.test(control$supportrelease_Rec3, notreat$supportrelease_Rec3)
tb.b <- t.test(t1$supportrelease_Rec3, notreat$supportrelease_Rec3)
tb.c <- t.test(t2$supportrelease_Rec3, notreat$supportrelease_Rec3) 

t.test(t1$supportrelease_Rec3, control$supportrelease_Rec3)
t.test(t2$supportrelease_Rec3, control$supportrelease_Rec3)
t.test(t2$supportrelease_Rec3, t1$supportrelease_Rec3)


##Recover adjusted p values using Holm correction method

pairwise.t.test(covid.white$supportrelease_Rec3, covid.white$condition, p.adjust.method="holm")

##Descriptives by treatment group
prop.table(table(notreat$supportrelease_Rec3))
prop.table(table(t1$supportrelease_Rec3))
prop.table(table(t2$supportrelease_Rec3))


treateffect.releasewrite <- c(0, tb.a$estimate[1]- tb.a$estimate[2], tb.b$estimate[1]- tb.b$estimate[2], tb.c$estimate[1]- tb.c$estimate[2])

lower.releasewrite <- c(0, tb.a$conf.int[1], tb.b$conf.int[1], tb.c$conf.int[1])

upper.releasewrite <- c(0, tb.a$conf.int[2], tb.b$conf.int[2], tb.c$conf.int[2])


treatment <- c("No\n treatment", "Informational\n Control",  "Perspective:\n egocentric", "Perspective:\n surrogate")

data.releasewrite <- as.data.frame(treateffect.releasewrite)
data.releasewrite$lower <- lower.releasewrite
data.releasewrite$upper <- upper.releasewrite
data.releasewrite$treatment <- treatment
data.releasewrite$treatment <- factor(data.releasewrite$treatment, levels = c("Perspective:\n surrogate", "Perspective:\n egocentric",  "Informational\n Control", "No\n treatment"))

data.releasewrite$labels <- formatC(data.releasewrite$treateffect.releasewrite, format = "f", digits = 3)
data.releasewrite$labels[data.releasewrite$labels == "0.000"] <- NA
data.releasewrite$lower[data.releasewrite$lower == 0] <- NA
data.releasewrite$upper[data.releasewrite$upper == 0] <- NA
data.releasewrite$labels2 <- c("", "0.060", "", "")


p2 <- ggplot(data = data.releasewrite, aes(x = treateffect.releasewrite, y = treatment)) +
  geom_text(aes(label = labels), nudge_y = 0.2, nudge_x = -0.002, size=6) +
  theme_bw() +
  geom_point(size = 4) +
  geom_errorbarh(aes(y = treatment, xmin = lower, xmax = upper, height=.1),size=.5) +
  xlab("Effect on Writing in Support")  + ylab("") +
  theme(axis.text = element_text(size = 18), #Increase size of axis text
        axis.title = element_text(size = 18))	+  geom_vline(xintercept = 0,size=.5,colour="black",linetype="dotted") +
  #ggtitle(title)+
  theme(legend.position = "none") +
  theme(plot.margin = unit(c(0.1, .5, 0.1, 0), "cm")) + 
  annotate("text", x = -0.07, y = 4.15, parse = TRUE, label = "bar(x)==-0.088", col = "black", size = 6) +
  xlim(-.2, .3)



pdf("fig2b.pdf", height= 5, width = 7)
p2
dev.off()



##################
##Who Responds?: Carceral Contact and Partisanship
##################


table(covid.white$pers_or_prox)
contact <- subset(covid.white, pers_or_prox == 1)
nocontact <- subset(covid.white, pers_or_prox == 0)
table(contact$pers_or_prox)
table(nocontact$pers_or_prox)

##means and effects
aggregate(nocontact$release, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE)

aggregate(nocontact$release, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(nocontact$release, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE)[1,2] ##diffs from notreat

aggregate(nocontact$supportrelease_Rec3, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE)

aggregate(nocontact$supportrelease_Rec3, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(nocontact$supportrelease_Rec, by = list(Category = nocontact$condition), FUN = mean, na.rm = TRUE)[1,2] ##diffs from notreat

##contact

aggregate(contact$release, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)

aggregate(contact$release, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(contact$release, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)[1,2] ##diffs from notreat

aggregate(contact$supportrelease_Rec3, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)

aggregate(contact$supportrelease_Rec3, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)[2] - aggregate(contact$supportrelease_Rec, by = list(Category = contact$condition), FUN = mean, na.rm = TRUE)[1,2] ##diffs from notreat




###OLS models to create plots


contact_att_lm <- lm(release ~ condition*pers_or_prox, data = covid.white)
summary(contact_att_lm)

contact_beh_lm <- lm(supportrelease_Rec3 ~ condition*pers_or_prox, data = covid.white)
summary(contact_beh_lm)


##Policy attitudes outcome

ef_contact_att <- effect("condition : pers_or_prox", contact_att_lm)

ef_contact_att

d_contact_att <- as.data.frame(ef_contact_att) #save as data frame
d_contact_att

#Save only estimates at 1 and 0
d_contact_att <- subset(d_contact_att, pers_or_prox == 0 | pers_or_prox == 1)
d_contact_att



#Add in factor names for plotting
d_contact_att$names <- recode(d_contact_att$pers_or_prox, '0' = 'No Contact' , '1' = 'Contact')
d_contact_att

d_contact_att$cond_lab <- recode(d_contact_att$condition, 'notreat' = 'No\nTreatment' , 'control' = 'Informational\nControl', 't1' = 'Perspective:\nEgocentric', 't2' = 'Perspective:\nSurrogate')
d_contact_att



############
##Figure 3a
############

contactplot <- ggplot(d_contact_att, aes(x = names, y = fit, shape = cond_lab, color = cond_lab)) +  
  scale_color_manual(values=c("#7db0ea", "#447fdd", "#225bb2", "#133e7e"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.3)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.3)) +
  theme_bw()  +
  ylab("Policy Opinion") + xlab(" ")+
  ylim(5, 12) + 
  scale_shape_manual(values = c(18, 15, 16, 17), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))

pdf("fig3a.pdf", height= 7, width = 8)
contactplot
dev.off()


## t tests for comparisons in text (policy attitudes): 

contact.nt <- subset(contact, condition == "notreat")
contact.c <- subset(contact, condition == "control")
contact.t1 <- subset(contact, condition == "t1")
contact.t2 <- subset(contact, condition == "t2")

t.test(contact.nt$release, contact.c$release)
t.test(contact.nt$release, contact.t1$release)
t.test(contact.nt$release, contact.t2$release)

t.test(contact.c$release, contact.t1$release)
t.test(contact.c$release, contact.t2$release)


nocontact.nt <- subset(nocontact, condition == "notreat")
nocontact.c <- subset(nocontact, condition == "control")
nocontact.t1 <- subset(nocontact, condition == "t1")
nocontact.t2 <- subset(nocontact, condition == "t2")

t.test(nocontact.nt$release, nocontact.c$release)
t.test(nocontact.nt$release, nocontact.t1$release)
t.test(nocontact.nt$release, nocontact.t2$release)

t.test(nocontact.c$release, nocontact.t1$release)
t.test(nocontact.c$release, nocontact.t2$release)




##Semibehavioral Outcome

ef_contact_beh <- effect("condition : pers_or_prox", contact_beh_lm)

ef_contact_beh

d_contact_beh <- as.data.frame(ef_contact_beh) #save as data frame
d_contact_beh

#Save only estimates at 1 and 0
d_contact_beh <- subset(d_contact_beh, pers_or_prox == 0 | pers_or_prox == 1)
d_contact_beh



#Add in factor names for plotting
d_contact_beh$names <- recode(d_contact_beh$pers_or_prox, '0' = 'No Contact' , '1' = 'Contact')
d_contact_beh

d_contact_beh$cond_lab <- recode(d_contact_beh$condition, 'notreat' = 'No\nTreatment' , 'control' = 'Informational\nControl', 't1' = 'Perspective:\nEgocentric', 't2' = 'Perspective:\nSurrogate')
d_contact_beh

############
##Figure 3b
############

contactplot_beh <- ggplot(d_contact_beh, aes(x = names, y = fit, shape = cond_lab, color = cond_lab)) +  
  scale_color_manual(values=c("#7db0ea", "#447fdd", "#225bb2", "#133e7e"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) + 
  theme_bw()  +
  ylab("Semibehavioral Outcome") + xlab(" ")+
  ylim(-.25, .25) + 
  scale_shape_manual(values = c(18, 15, 16, 17), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("fig3b.pdf", height= 7, width = 8)
contactplot_beh
dev.off()



## t tests for comparisons in text (semibehavioral): 

t.test(contact.nt$supportrelease_Rec3, contact.c$supportrelease_Rec3)
t.test(contact.nt$supportrelease_Rec3, contact.t1$supportrelease_Rec3)
t.test(contact.nt$supportrelease_Rec3, contact.t2$supportrelease_Rec3)

t.test(contact.c$supportrelease_Rec3, contact.t1$supportrelease_Rec3)
t.test(contact.c$supportrelease_Rec3, contact.t2$supportrelease_Rec3)

t.test(nocontact.nt$supportrelease_Rec3, nocontact.c$supportrelease_Rec3)
t.test(nocontact.nt$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(nocontact.nt$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)

t.test(nocontact.c$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(nocontact.c$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)



##comparing across contact/nocontact
t.test(contact.nt$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(contact.nt$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)

t.test(contact.c$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(contact.c$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)

t.test(contact.t1$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)
t.test(contact.t1$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)
t.test(contact.t2$supportrelease_Rec3, nocontact.t2$supportrelease_Rec3)
t.test(contact.t2$supportrelease_Rec3, nocontact.t1$supportrelease_Rec3)



#########
## Partisanship 
#########


table(covid.white$pid7)
dems <- subset(covid.white, pid7 < 3)
table(dems$pid7)
reps <- subset(covid.white, pid7 > 3)
table(reps$pid7)

aggregate(dems$release, by = list(Category = dems$condition), FUN = mean, na.rm = TRUE)
aggregate(reps$release, by = list(Category = reps$condition), FUN = mean, na.rm = TRUE) 


covid.white$dem[covid.white$pid7 < 3] <- 1
covid.white$dem[covid.white$pid7 > 3] <- 0
table(covid.white$pid7, covid.white$dem)


##OLS models for creating figure 4
pid_att_lm <- lm(release ~ condition*dem, data = covid.white)
summary(pid_att_lm)


pid_beh_lm <- lm(supportrelease_Rec3 ~ condition*dem, data = covid.white)
summary(pid_beh_lm)


ef_pid_att <- effect("condition : dem", pid_att_lm)

ef_pid_att

d_pid_att <- as.data.frame(ef_pid_att) #save as data frame
d_pid_att

#Save only estimates at 1 and 0
d_pid_att <- subset(d_pid_att, dem == 0 | dem == 1)
d_pid_att



#Add in factor names for plotting
d_pid_att$names <- recode(d_pid_att$dem, '0' = 'Republicans' , '1' = 'Democrats')
d_pid_att

d_pid_att$cond_lab <- recode(d_pid_att$condition, 'notreat' = 'No\nTreatment' , 'control' = 'Informational\nControl', 't1' = 'Perspective:\nEgocentric', 't2' = 'Perspective:\nSurrogate')
d_pid_att

###########
##Figure 4a
###########

pidplot <- ggplot(d_pid_att, aes(x = names, y = fit, shape = cond_lab, color = cond_lab)) +  
  scale_color_manual(values=c("#7db0ea", "#447fdd", "#225bb2", "#133e7e"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) +
  theme_bw()  +
  ylab("Policy Opinion") + xlab(" ")+
  ylim(5, 13) + 
  scale_shape_manual(values = c(18, 15, 16, 17), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("fig4a.pdf", height= 7, width = 8)
pidplot
dev.off()


##t tests for comparisons in text (policy attitudes)

dems.nt <- subset(dems, condition == "notreat")
dems.c <- subset(dems, condition == "control")
dems.t1 <- subset(dems, condition == "t1")
dems.t2 <- subset(dems, condition == "t2")

t.test(dems.nt$release, dems.c$release)
t.test(dems.nt$release, dems.t1$release)
t.test(dems.nt$release, dems.t2$release)

t.test(dems.c$release, dems.t1$release)
t.test(dems.c$release, dems.t2$release)


reps.nt <- subset(reps, condition == "notreat")
reps.c <- subset(reps, condition == "control")
reps.t1 <- subset(reps, condition == "t1")
reps.t2 <- subset(reps, condition == "t2")

t.test(reps.nt$release, reps.c$release)
t.test(reps.nt$release, reps.t1$release)
t.test(reps.nt$release, reps.t2$release)

t.test(reps.c$release, reps.t1$release)
t.test(reps.c$release, reps.t2$release)




### Semibehavioral outcome

aggregate(dems$supportrelease_Rec3, by = list(Category = dems$condition), FUN = mean, na.rm = TRUE)
aggregate(reps$supportrelease_Rec3, by = list(Category = reps$condition), FUN = mean, na.rm = TRUE) 

ef_pid_beh <- effect("condition : dem", pid_beh_lm)

ef_pid_beh

d_pid_beh <- as.data.frame(ef_pid_beh) #save as data frame
d_pid_beh

#Save only estimates at 1 and 0
d_pid_beh <- subset(d_pid_beh, dem == 0 | dem == 1)
d_pid_beh



#Add in factor names for plotting
d_pid_beh$names <- recode(d_pid_beh$dem, '0' = 'Republicans' , '1' = 'Democrats')
d_pid_beh

d_pid_beh$cond_lab <- recode(d_pid_beh$condition, 'notreat' = 'No\nTreatment' , 'control' = 'Informational\nControl', 't1' = 'Perspective:\nEgocentric', 't2' = 'Perspective:\nSurrogate')
d_pid_beh


############
##Figure 4b
############

pidplot_beh <- ggplot(d_pid_beh, aes(x = names, y = fit, shape = cond_lab, color = cond_lab)) +  
  scale_color_manual(values=c("#7db0ea", "#447fdd", "#225bb2", "#133e7e"), name=" ") + 
  geom_point(size = 4, position = position_dodge(.2)) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.15, position = position_dodge(.2)) + 
  theme_bw()  +
  ylab("Semibehavioral Outcome") + xlab(" ")+
  ylim(-.3, .3) + 
  scale_shape_manual(values = c(18, 15, 16, 17), name=" ")+
  theme(legend.position = "bottom", legend.text = element_text(size = 16),
        axis.text.x = element_text(size = 16), 
        axis.text.y = element_text(size = 18), 
        axis.title.y = element_text(size = 18))


pdf("fig4b.pdf", height= 7, width = 8)
pidplot_beh
dev.off()





##t tests for comparisons in text (Semibehavioral)


t.test(dems.nt$supportrelease_Rec3, dems.c$supportrelease_Rec3)
t.test(dems.nt$supportrelease_Rec3, dems.t1$supportrelease_Rec3)
t.test(dems.nt$supportrelease_Rec3, dems.t2$supportrelease_Rec3)

t.test(dems.c$supportrelease_Rec3, dems.t1$supportrelease_Rec3)
t.test(dems.c$supportrelease_Rec3, dems.t2$supportrelease_Rec3)


t.test(reps.nt$supportrelease_Rec3, reps.c$supportrelease_Rec3)
t.test(reps.nt$supportrelease_Rec3, reps.t1$supportrelease_Rec3)
t.test(reps.nt$supportrelease_Rec3, reps.t2$supportrelease_Rec3)

t.test(reps.c$supportrelease_Rec3, reps.t1$supportrelease_Rec3)
t.test(reps.c$supportrelease_Rec3, reps.t2$supportrelease_Rec3)


t.test(dems.nt$supportrelease_Rec3, reps.t1$supportrelease_Rec3)
t.test(dems.nt$supportrelease_Rec3, reps.t2$supportrelease_Rec3)


